/* (c) Copyright 1993 by Joseph Felsenstein.
   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */

#include <math.h>
#include <map>
#include <iostream>
#include <sstream>
#include <algorithm>

// component headers
#include "tree_ops.h"
#include "matrix.h"
#include "sparse_matrix.h"
#include "tree_kitsch.h"
#include "tree_kitsch.priv.h"
#include "tree_maker_nh.h"
#include "tree_maker_tree.h"
#include "digest.h"
#include "md4.h"
#define DEBUG_LEVEL 2
#include "tree_UPGMA.h"
#include "debug_func.h"

using std::vector;
using std::string;
using std::copy;
using std::cerr;
using std::cout;
using std::ostringstream;
using bcpp_tree::t_tree;
using bcpp_tree::t_node;



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
// globals
namespace
{
	inline	float square(const float d){return d * d;}
}
//
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



//	for (unsigned i = 1; i < cnt_leaves; ++i)
//		for (unsigned j = 0; j < i; ++j)
//	for (unsigned i = 0; i < cnt_leaves - 1 ; ++i)
//		for (unsigned j = i + 1; j < cnt_leaves; ++j)
//	two forms are identical





namespace kitsch_tree
{


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	struct t_kitsch_node

//

#define epsilon         0.000001   /* a very small but not too small number */


//________________________________________________________________________________________

//	CHECK_LINKS

//________________________________________________________________________________________
void t_kitsch_node::CHECK_LINKS() const
{
#ifndef NDEBUG
	static unsigned indent = 0;
	if (!indent)
		debug_cout(5, "CHECK_LINKS node ENTRY");
	debug_cout(5, string(indent, ' ') << get_UID());
	indent += 2;
	if (get_parent())
		assert(get_parent()->children[child_num] == this);
	if (!is_leaf())
	{
		if (children[CHILD_LEFT])
		{
			assert(children[CHILD_LEFT]->get_parent() == this);
			assert(children[CHILD_LEFT]->child_num == CHILD_LEFT);
			children[CHILD_LEFT]->CHECK_LINKS();
		}
		if (children[CHILD_RITE])
		{
			assert(children[CHILD_RITE]->get_parent() == this);
			assert(children[CHILD_RITE]->child_num == CHILD_RITE);
			children[CHILD_RITE]->CHECK_LINKS();
		}
	}
	indent -= 2;
	if (!indent)
		debug_cout(5, "CHECK_LINKS node EXIT");
#endif
}


//________________________________________________________________________________________

//	sibling

//________________________________________________________________________________________
inline t_kitsch_node* t_kitsch_node::sibling() const
{
	assert (get_parent());
	return get_parent()->children[(child_num + 1) % CHILD__ALL];
}
//________________________________________________________________________________________

//	set_parent

//________________________________________________________________________________________
inline void t_kitsch_node::set_parent(t_kitsch_node* new_parent, eCHILD_NUM num)
{
#ifndef NDEBUG
	debug_func_cout(5);
#endif

	// old parent's children should be null
	if (get_parent())
	{
		assert(get_parent()->children[child_num] == this);
		get_parent()->children[child_num] = 0;
	}

	parent = new_parent;

	if (new_parent)
	{
		// new parent's old children should be orphaned
		if (new_parent->children[num])
			new_parent->children[num]->parent = 0;
		new_parent->children[num] = this;
	}
	child_num = num;
}


//________________________________________________________________________________________

//	assignment operator =

//________________________________________________________________________________________
t_kitsch_node& t_kitsch_node::operator = (const t_kitsch_node& other)
{
	child_num			= other.child_num;
	len_to_leaf			= other.len_to_leaf;
	same_len_to_leaf	= other.same_len_to_leaf;
	weight				= other.weight;
	processed			= other.processed;
	UID					= other.get_UID();
//	treepos				= other.treepos;
	leaf				= other.leaf;
	return *this;
}

//________________________________________________________________________________________

//	operator==

//________________________________________________________________________________________
bool t_kitsch_node::operator!= (const t_kitsch_node& other) const
{
	return !(*this == other);
}

bool t_kitsch_node::operator== (const t_kitsch_node& other) const
{
	return (
			len_to_leaf			== other.len_to_leaf		&&
			same_len_to_leaf	== other.same_len_to_leaf	&&
			weight				== other.weight				&&
			processed			== other.processed			&&
			UID					== other.get_UID()				&&
			treepos				== other.treepos			);
}


//________________________________________________________________________________________

//	left_splice_below

//	inserts the nodes new_parent and its left descendant, newtip to the tree.
//	at becomes new_parent's right descendant

//________________________________________________________________________________________
inline void t_kitsch_node::left_splice_below(t_kitsch_node& new_parent, t_kitsch_node& new_tip)
{
#ifndef NDEBUG
	debug_func_cout(5);
	if (get_parent())
		debug_cout(5, get_parent()->get_UID() << " = new parent");
	debug_cout(5, new_parent.get_UID() << " = new top");
	debug_cout(5, new_tip.get_UID() << " left");
	debug_cout(5, get_UID() << " right");
	tree->CHECK_LINKS();
#endif



	// make new_parent children of parent in my place
	new_parent.set_parent(get_parent(), child_num);

	// now link up this and new_tip as children of new_parent
	new_tip.set_parent(&new_parent, CHILD_LEFT);
	set_parent(&new_parent, CHILD_RITE);


	// reroot if necessary
	if (tree->root == this)
		tree->root = &new_parent;
	tree->root->parent = 0;




#ifndef NDEBUG
	tree->CHECK_LINKS();
#endif


}  /* left_splice_below */



//

//	struct t_kitsch_node

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888








//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


//	operations
//	=============


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//________________________________________________________________________________________

//	op_print_tree


//________________________________________________________________________________________
class op_print_tree
{
	std::ostream&		os;
	t_kitsch_tree&		tree;
	unsigned			precision;
	unsigned			curr_level;
	std::vector<char>	levels;
public:
	op_print_tree(t_kitsch_tree& tree_, unsigned p, std::ostream& os_):
			 os(os_), tree(tree_), precision(p), curr_level(0){}
	void print_horizontal_lines()
	{
		os << std::string(precision + 1, '-');
	}
	void print_vertical_lines()
	{
		for (unsigned i = 0; i < levels.size(); ++i)
			os << std::string(precision, ' ') << (levels[i] ? '|' : ' ');
	}
	void print(t_kitsch_node& curr_node)
	{
		recurse(curr_node);
	}
	void recurse(t_kitsch_node& curr_node);
};

void op_print_tree::recurse(t_kitsch_node& curr_node)
{
	if (curr_node.is_leaf())
	{
		print_horizontal_lines();
		os << "[" << curr_node.get_UID() << "] " << tree.names[curr_node.get_UID()] << std::endl;
		return;
	}
	else
	{
		print_horizontal_lines();
		levels.push_back(1);
		recurse(*curr_node.children[CHILD_LEFT]);

		print_vertical_lines();
 		os << "[" << curr_node.get_UID() << "] " << std::endl;

		levels.back() = 0;

		print_vertical_lines();
		recurse(*curr_node.children[CHILD_RITE]);

		levels.pop_back();
	}
}







//________________________________________________________________________________________

//	op_collect_nodes

//		go through a tree and collect all children nodes

//________________________________________________________________________________________
class op_collect_nodes
{
	vector<t_kitsch_node*> leaf_nodes;
public:
	vector<t_kitsch_node*> nodes;
	unsigned& cnt_leaves;
	unsigned& cnt_nodes;
	op_collect_nodes(unsigned& s, unsigned& n): cnt_leaves(s), cnt_nodes(n){}
	void operator()(t_kitsch_node& curr_node)
	{
		nodes.clear();
		recurse(curr_node);
		cnt_leaves = leaf_nodes.size();
		// leaf nodes come first before internal nodes
		leaf_nodes.insert(leaf_nodes.end(), nodes.begin(), nodes.end());
		leaf_nodes.swap(nodes);
		cnt_nodes = nodes.size();
		assert(cnt_leaves * 2 - 1 == cnt_nodes);
	}
	void recurse(t_kitsch_node& curr_node)
	{
		if (curr_node.is_leaf())
			leaf_nodes.push_back(&curr_node);
		else
		{
			nodes.push_back(&curr_node);
			recurse(*curr_node.children[CHILD_LEFT]);
			recurse(*curr_node.children[CHILD_RITE]);
		}
	}
};

//________________________________________________________________________________________

//	op_optimize_sub_trees

//		go through a tree and optimize all branches with fewer leaves than max_tree_size

//________________________________________________________________________________________
class op_optimize_sub_trees
{
    t_kitsch_tree&	tree;
	unsigned		njumble;
	unsigned		max_tree_size;
	unsigned&		cnt_trees_examined;
public:
	op_optimize_sub_trees(	t_kitsch_tree& tree_,
							unsigned njumble_,
							unsigned max_tree_size_,
							unsigned& cnt_trees_examined_):
			 tree(tree_), njumble(njumble_), max_tree_size(max_tree_size_),
			cnt_trees_examined(cnt_trees_examined_){}

	// use pointer because it might change
	void optimize(t_kitsch_node& curr_node);
	void start_at(t_kitsch_node& curr_node);
	bool recurse(t_kitsch_node& curr_node, unsigned& branch_leaf_count);
};

void op_optimize_sub_trees::optimize(t_kitsch_node& curr_node)
{
//	op_print_tree(tree, 6, cout).print(*tree.root);
	t_kitsch_tree best_tree(tree, curr_node);
	//cerr << "\n" << best_tree.cnt_leaves << "\n";
	vector<unsigned> order_of_entry;
	for (unsigned i = 0; i < best_tree.cnt_leaves; ++i)
		order_of_entry.push_back(i);
	for (unsigned i = 0; i < njumble; ++i)
	{
		t_kitsch_tree curr_tree(tree, curr_node);
		std::reverse(order_of_entry.begin(), order_of_entry.end());
		for (unsigned j = 0; j < order_of_entry.size(); ++j)
			std::swap(order_of_entry[j], order_of_entry[rand() % order_of_entry.size()]);

		curr_tree.maketree(order_of_entry);
		if (curr_tree.best_sum_sqrs > best_tree.best_sum_sqrs)
			best_tree = curr_tree;
		cnt_trees_examined += curr_tree.cnt_trees_examined;
	#ifndef NDEBUG
			best_tree.CHECK_LINKS();
			curr_tree.CHECK_LINKS();
	#endif

	}
//	op_print_tree(best_tree, 6, cout).print(*best_tree.root);
	tree.replace_subtree(best_tree, curr_node);
	tree.CHECK_LINKS();
	tree.pdots->print_indicator('!');

//	op_print_tree(tree, 6, cout).print(*tree.root);
//	cout << string (60, '8') << std::endl << std::endl;
}

void op_optimize_sub_trees::start_at(t_kitsch_node& curr_node)
{
	unsigned leaf_count = 0;
	if (!recurse(curr_node, leaf_count))
		optimize(curr_node);
	tree.evaluate();
}

bool op_optimize_sub_trees::recurse(t_kitsch_node& curr_node, unsigned& branch_leaf_count)
{
	if (curr_node.is_leaf())
	{
		++branch_leaf_count;
		return false;
	}

	// if internal node, get count of all descendent leaves
	//		and add that to the parent's running total
	unsigned children_leaf_count_l = 0;
	unsigned children_leaf_count_r = 0;
	bool left_optimized = recurse(*curr_node.children[CHILD_LEFT], children_leaf_count_l);
	bool rite_optimized = recurse(*curr_node.children[CHILD_RITE], children_leaf_count_r);
	branch_leaf_count += children_leaf_count_l;
	branch_leaf_count += children_leaf_count_r;

	// optimize the other branch if only one branch optimized
	if (left_optimized)
	{
		if (!rite_optimized && children_leaf_count_r > 4)
		{
			//cerr << "\n" << children_leaf_count_r << "\n";
			optimize(*curr_node.children[CHILD_RITE]);
		}
		return true;
	}
	else if (rite_optimized)
	{
		if (children_leaf_count_l > 4)
		{
			//cerr << "\n" << children_leaf_count_l << "\n";
			optimize(*curr_node.children[CHILD_LEFT]);
		}
		return true;
	}

	// if reach threshold, optimize the branch starting with me
	else if (children_leaf_count_l + children_leaf_count_r >= max_tree_size)
	{
		if (children_leaf_count_l <=5 || children_leaf_count_r <=5)
		{
			//cerr << "\n" << children_leaf_count_l << "+" << children_leaf_count_r << " b\n";
			optimize(curr_node);
		}
		else
		{
			//cerr << "\n" << children_leaf_count_l << "+" << children_leaf_count_r << "\n";
			optimize(*curr_node.children[CHILD_LEFT]);
			optimize(*curr_node.children[CHILD_RITE]);
		}
		return true;
	}

	// nothing optimized
	return false;
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


//	operations
//	=============


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888













//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


//	t_kitsch_tree
//	=============


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//________________________________________________________________________________________

//	c'tor

//________________________________________________________________________________________
t_kitsch_tree::t_kitsch_tree(	bcpp_matrix::t_dense_matrix_f&	d_,
								bcpp_matrix::t_dense_matrix_f&	w_,
								const vector<string >&  names_,
								float					power_,
								t_progress_indicator& 	dots_,
								unsigned				cnt_leaves_)
	:d(d_), w(w_), names(names_), power(power_), cnt_leaves(cnt_leaves_),
	cnt_trees_examined(0),
	cnt_non_null_data(0), sum_sqrs(-1e100), best_sum_sqrs(-1e100), pdots(&dots_),
	root(0), cached_evaluate_results(10000), cnt_evaluate_cache_hits(0)
{
#ifndef NDEBUG
	debug_func_cout(5);
#endif
	if (!cnt_leaves)
		cnt_leaves = names.size();
	cnt_nodes = cnt_leaves * 2 - 1;
	assert(d.size() >= cnt_nodes);
	assert(w.size() >= cnt_nodes);

	// allocate memory for actual nodes
	nulls.allocate(cnt_leaves, false);
	shared_mem.assign(cnt_leaves,
							t_node_type(cnt_nodes, this, t_kitsch_node::eLEAF));
	shared_mem.insert(shared_mem.end(),  (cnt_leaves - 1),
							t_node_type(cnt_nodes, this, t_kitsch_node::eNON_LEAF));

	// all pointers set to zero
	nodes.assign(cnt_nodes, 0);

	// all pointers set to allocated memory and set UID of nodes
	// indices are immutable
	for (unsigned i = 0; i < nodes.size(); ++i)
	{
		nodes[i] = &shared_mem[i];
		nodes[i]->UID		=
		nodes[i]->treepos	= i;
	}

#if (!defined(NDEBUG))
	// mark leaves
	for (unsigned i = 0; i < cnt_leaves; ++i)
		assert(nodes[i]->is_leaf());
	// mark leaves
	for (unsigned i = cnt_leaves; i < shared_mem.size(); ++i)
		assert(!nodes[i]->is_leaf());
#endif
}


//________________________________________________________________________________________

//	c'tor 2

//		create tree from a specified sub-branch of another tree

//________________________________________________________________________________________
t_kitsch_tree::t_kitsch_tree(t_kitsch_tree other, t_node_type& from_node)
	:d(other.d), w(other.w), names(other.names), power(other.power),
	cnt_trees_examined(0),
	cnt_non_null_data(0), sum_sqrs(-1e100), best_sum_sqrs(-1e100), pdots(other.pdots),
	root(0), cnt_evaluate_cache_hits(0)
{
#ifndef NDEBUG
	debug_func_cout(5);
#endif

	// collect all descendent nodes with leaf nodes before internal nodes
	op_collect_nodes subtree(cnt_leaves, cnt_nodes);
	subtree(from_node);

	assert(d.size() >= cnt_nodes);
	assert(w.size() >= cnt_nodes);

	// allocate memory for actual nodes
	nulls.allocate(cnt_leaves, false);
	shared_mem.assign(cnt_leaves,
							t_node_type(cnt_nodes, this, t_kitsch_node::eLEAF));
	shared_mem.insert(shared_mem.end(),  (cnt_leaves - 1),
							t_node_type(cnt_nodes, this, t_kitsch_node::eNON_LEAF));

	// all pointers set to zero
	nodes.assign(cnt_nodes, 0);

	// all pointers set to allocated memory and set UID of nodes
	// indices are immutable
	for (unsigned i = 0; i < cnt_nodes; ++i)
	{
		nodes[i] = &shared_mem[i];
		// UID is the same UID as the subtree
		*nodes[i]			= *subtree.nodes[i];
		// but node is my own
		nodes[i]->treepos	= i;
	}

#if (!defined(NDEBUG))
	// mark leaves
	for (unsigned i = 0; i < cnt_leaves; ++i)
		assert(nodes[i]->is_leaf());
	// mark leaves
	for (unsigned i = cnt_leaves; i < shared_mem.size(); ++i)
		assert(!nodes[i]->is_leaf());
#endif

}


//________________________________________________________________________________________

//	replace_subtree

//		replace a tree branch with an entire other tree

//		this is intended to replace one branch with a tree with the same number of nodes
//		i.e. the total number of nodes etc. should not change!!!

//________________________________________________________________________________________
void t_kitsch_tree::replace_subtree(t_kitsch_tree other, t_node_type& from_node)
{
#ifndef NDEBUG
	debug_func_cout(5);
#endif

	// save from_node parent
	t_node_type* from_node_parent = from_node.get_parent();
	eCHILD_NUM from_node_child_num = from_node.child_num;
	if (from_node_parent)
		from_node_parent->children[from_node_child_num] = 0;
	bool set_new_root = root == &from_node;

	// collect all descendent nodes with leaf nodes before internal nodes
	unsigned subtree_cnt_leaves, subtree_cnt_nodes;
	op_collect_nodes subtree(subtree_cnt_leaves, subtree_cnt_nodes);
	subtree(from_node);

	// only replace subtree with a tree of equivalent size
	assert(other.cnt_leaves == subtree_cnt_leaves);
	assert(other.cnt_nodes = subtree_cnt_nodes);


	// copy nodes
	for (unsigned i = 0; i < other.cnt_nodes; ++i)
		*subtree.nodes[i] = *other.nodes[i];


	for (unsigned i = 0; i < other.cnt_nodes; ++i)
	{
		subtree.nodes[i]->parent = 0;
		subtree.nodes[i]->children[CHILD_LEFT] =
		subtree.nodes[i]->children[CHILD_RITE] = 0;
	}

		// restore links to parent
	for (unsigned i = 0; i < other.nodes.size(); ++i)
	{
		if (other.nodes[i] == other.root)
		{
			subtree.nodes[i]->set_parent(from_node_parent, from_node_child_num);
			if (set_new_root)
				root = subtree.nodes[i];
			continue;
		}

		assert (other.nodes[i]->parent->treepos < subtree.cnt_nodes);
		subtree.nodes[i]->set_parent(subtree.nodes[other.nodes[i]->get_parent()->treepos],
										subtree.nodes[i]->child_num);
	}
}

//________________________________________________________________________________________

//	prepare_data

//		assign data to weights and prepare distances
//		only need to do once for all trees which share the same matrix and weights

//________________________________________________________________________________________
void t_kitsch_tree::prepare_data()
{
#ifndef NDEBUG
	debug_func_cout(5);
	assert(d.size() >= cnt_nodes);
	assert(w.size() >= cnt_nodes);
	assert(cnt_nodes == nodes.size());
	assert(cnt_leaves* 2 - 1 == cnt_nodes);
#endif

	// assign data to nodes distances and weights
	if (power == 2.0)
	{
		for (unsigned i = 0; i < cnt_leaves - 1 ; ++i)
			for (unsigned j = i + 1; j < cnt_leaves; ++j)
			{
				unsigned ii = nodes[i]->get_UID();
				unsigned jj = nodes[j]->get_UID();
				// nulls
				if (d[ii][jj] > 5000.0 || d[ii][jj] < 0.0)
				{
					nulls[ii][jj] = nulls[jj][ii] = true;
					w[ii][jj] = w[jj][ii] = 0.0;
				}
				else
				{
					if (d[ii][jj] < epsilon)
						d[ii][jj] =	d[jj][ii] = epsilon;
					w[ii][jj] = w[jj][ii] =	1.0 / square (d[ii][jj]);
					++cnt_non_null_data;
				}
			}
	}
	else
	{
		for (unsigned i = 0; i < cnt_leaves - 1 ; ++i)
			for (unsigned j = i + 1; j < cnt_leaves; ++j)
			{
				unsigned ii = nodes[i]->get_UID();
				unsigned jj = nodes[j]->get_UID();
				// nulls
				if (d[ii][jj] > 5000.0 || d[ii][jj] < 0.0)
				{
					nulls[ii][jj] = nulls[jj][ii] = true;
					w[ii][jj] = w[jj][ii] = 0.0;
				}
				else
				{
					if (d[ii][jj] < epsilon)
						d[ii][jj] =	d[jj][ii] = epsilon;
					w[ii][jj] = w[jj][ii] =	w[jj][ii] =	1.0 / pow (d[ii][jj], power);
					++cnt_non_null_data;
				}
			}

	}
}








//________________________________________________________________________________________

//	CHECK_LINKS()

//		for debugging

//________________________________________________________________________________________
void t_kitsch_tree::CHECK_LINKS() const
{
#ifndef NDEBUG
	debug_cout(5, "CHECK_LINKS ENTRY");
	for (unsigned i = 0; i < nodes.size(); ++i)
	{
		t_node_type* a = nodes[i];
		assert(a);
		if (a->get_parent() && a->get_parent()->children[a->child_num] != a)
		{
			cout << debug_prefix << a->get_UID()  << " failed " << (a->child_num == CHILD_LEFT ?
						 "left\n" : "right") << std::endl;
			cout << debug_prefix  << a->get_parent()->get_UID()  << " parent"  << std::endl;
			if (a->get_parent()->children[a->child_num])
				cout << debug_prefix  << a->get_parent()->children[a->child_num]->get_UID() << " children" << std::endl;
			assert(a->get_parent()->children[a->child_num] == a);
		}
		if (a->children[CHILD_LEFT])
		{
			if (a->children[CHILD_LEFT]->get_parent() != a)
			{
				cout << debug_prefix  << a->get_UID()  << " failed" << std::endl;
				cout << debug_prefix << a->children[CHILD_LEFT]->get_UID() << std::endl;
				if (a->children[CHILD_LEFT]->get_parent())
					cout << debug_prefix << a->children[CHILD_LEFT]->get_UID() << std::endl;
				assert(a->children[CHILD_LEFT]->get_parent() == a);
			}
			if (a->children[CHILD_LEFT]->child_num != CHILD_LEFT)
			{
				cout << debug_prefix << a->get_UID()  << " failed" << std::endl;
				cout << debug_prefix << a->children[CHILD_LEFT]->get_UID() << std::endl;
				cout << debug_prefix << (a->children[CHILD_LEFT]->child_num == CHILD_LEFT ?
						 "left\n" : "right") << " should be CHILD_LEFT" << std::endl;
				assert(a->children[CHILD_LEFT]->child_num == CHILD_LEFT);
			}
		}
		if (a->children[CHILD_RITE])
		{
			if (a->children[CHILD_RITE]->get_parent() != a)
			{
				cout << debug_prefix << a->get_UID()  << " failed" << std::endl;
				cout << debug_prefix << a->children[CHILD_RITE]->get_UID() << std::endl;
				if (a->children[CHILD_RITE]->get_parent())
					cout << debug_prefix << a->children[CHILD_RITE]->get_UID() << std::endl;
				assert(a->children[CHILD_RITE]->get_parent() == a);
			}
			if (a->children[CHILD_RITE]->child_num != CHILD_RITE)
			{
				cout << debug_prefix << a->get_UID()  << " failed" << std::endl;
				cout << debug_prefix << a->children[CHILD_RITE]->get_UID() << std::endl;
				cout << debug_prefix << (a->children[CHILD_RITE]->child_num == CHILD_RITE ?
						 "RITE\n" : "right") << "should be CHILD_RITE" << std::endl;
				assert(a->children[CHILD_RITE]->child_num == CHILD_RITE);
			}
		}

	}
	if (root)
		root->CHECK_LINKS();
	debug_cout(5, "CHECK_LINKS EXIT");
#endif
}



//________________________________________________________________________________________

//	remove

//	removes nodes item and its ancestor, parent, from the tree.
//		the new descendant of parent's parent (in the place of parent)
//		will be item's sibling
//		returns pointers to the deleted nodes, item and parent

//________________________________________________________________________________________
inline t_kitsch_tree::t_node_type* t_kitsch_tree::remove(t_node_type& item)
{
#ifndef NDEBUG
	CHECK_LINKS();
	debug_func_cout(5);
	debug_cout(5, item.get_UID() << " remove item");
	debug_cout(5, item.get_parent()->get_UID() << " and remove parent");
	if (item.get_parent()->get_parent())
		debug_cout(5, item.get_parent()->get_parent()->get_UID() << " from grandparent");
	if (item.sibling())
		debug_cout(5, item.sibling()->get_UID() << " sibling will be in the place of item");
#endif


	// forget it if has no parent: e.g. root
	t_node_type* parent = item.get_parent();
	if (!parent)
		return parent;

	t_kitsch_node* sibling = item.sibling();



	// reroot if necessary to the other children of root
	if (root == parent)
		root = sibling;
	if (sibling)
		sibling->set_parent(parent->get_parent(), parent->child_num);

	// break up links between item and parent as well
	item.parent		=0;
	parent->children[0]	=
	parent->children[1]	=0;



#ifndef NDEBUG
	assert(!parent->get_parent());
	root->CHECK_LINKS();
#endif
	return parent;
}  /* remove */




//________________________________________________________________________________________

//	scrunchtraverse

//		traverse to find closest t_kitsch_node to the current one

//________________________________________________________________________________________
void scrunchtraverse(t_kitsch_node *u, t_kitsch_node*& closest, float& len_to_leaf_max)
{
#ifndef NDEBUG
	debug_func_cout(4);
	assert(u);
#endif


	if (!u->same_len_to_leaf)
	{
		if (u->len_to_leaf > len_to_leaf_max)
		{
			closest = u;
			len_to_leaf_max = u->len_to_leaf;
		}
		return;
	}
	if (u->get_parent())
		u->len_to_leaf = u->get_parent()->len_to_leaf;
	if (!u->is_leaf())
	{
		scrunchtraverse(u->children[CHILD_LEFT], closest, len_to_leaf_max);
		scrunchtraverse(u->children[CHILD_RITE], closest, len_to_leaf_max);
	}
#ifndef NDEBUG
	debug_cout(5, "item = " << u->get_UID() << "  len_to_leaf_max=" << len_to_leaf_max);
#endif
}  /* scrunchtraverse */




//________________________________________________________________________________________

//	combine

//________________________________________________________________________________________
inline void combine(t_kitsch_node& a, t_kitsch_node& b)
{
#ifndef NDEBUG
	debug_func_cout(4);
#endif


	/* put t_kitsch_node b into the set having the same len_to_leaf as a */
	if (a.weight + b.weight <= 0.0)
		a.len_to_leaf = 0.0;
	else
		a.len_to_leaf = (a.len_to_leaf * a.weight + b.len_to_leaf * b.weight) /
						(a.weight + b.weight);
	a.weight += b.weight;
	b.same_len_to_leaf = true;
}  /* combine */




//________________________________________________________________________________________

//	scrunch

//		see if nodes can be combined to prevent negative lengths

//________________________________________________________________________________________
void scrunch(t_kitsch_node& s)
{
#ifndef NDEBUG
	debug_func_cout(3);
#endif
	t_kitsch_node *closest	= 0;
	do
	{
		float len_to_leaf_max	= -1.0;
		if (!s.is_leaf())
		{
			scrunchtraverse(s.children[CHILD_LEFT], closest, len_to_leaf_max);
			scrunchtraverse(s.children[CHILD_RITE], closest, len_to_leaf_max);
		}
		if (len_to_leaf_max <= s.len_to_leaf)
			break;
		assert(closest);
		combine(s, *closest);
	} while (1);
}  /* scrunch */




//________________________________________________________________________________________

//	secondtraverse

//		recalculate distances, add to sum

//________________________________________________________________________________________
float t_kitsch_tree::secondtraverse(	t_kitsch_node& ll,
										t_kitsch_node& kk,
										t_kitsch_node& ii,
										t_kitsch_node& jj)
{
#ifndef NDEBUG
	debug_func_cout(4);
	debug_cout(5, ll.get_UID());
#endif

	float sum = 0.0;
	if (!(ll.processed || ll.is_leaf()))
	{
		sum += secondtraverse(*ll.children[CHILD_LEFT], kk,ii,jj);
		sum += secondtraverse(*ll.children[CHILD_RITE], kk,ii,jj);
		return sum;
	}
	if (!(ll.get_UID() != kk.get_UID() && ll.processed))
		return sum;


	unsigned i = ii.get_UID();
	unsigned j = jj.get_UID();
	unsigned k = kk.get_UID();
	unsigned l = ll.get_UID();
	assert(i < w.size());
	assert(j < w.size());
	assert(k < w.size());
	assert(l < w.size());

// symmetrical matrices
	float wil = w[i][l];
	float wjl = w[j][l];
	float wkl = wil + wjl;
	w[k][l] = w[l][k] = wkl;

	if (wkl <= 0.0)
		d[k][l] =
		d[l][k] = 0.0;
	else
		d[k][l]	=
		d[l][k]	= (wil * d[i][l] + wjl * d[j][l]) / wkl;
	if (wkl > 0.0)
		sum += 2 * wil * wjl / wkl * square(d[i][l] - d[j][l]);
	return sum;
}  /* secondtraverse */




//________________________________________________________________________________________

//	firstraverse

//		go through tree calculating branch lengths

//________________________________________________________________________________________
float t_kitsch_tree::firstraverse(t_kitsch_node& q, t_kitsch_node& r)
{
#ifndef NDEBUG
	debug_func_cout(3);
	debug_cout(5, q.get_UID());
	assert(q.is_leaf() || q.children[CHILD_LEFT]);
	if (!q.is_leaf())
		debug_cout(5, q.children[CHILD_LEFT]->get_UID() << " left");
	assert(q.is_leaf() || q.children[CHILD_RITE]);
	if (!q.is_leaf())
		debug_cout(5, q.children[CHILD_RITE]->get_UID() << " right");
#endif


	float sum = 0.0;

	q.same_len_to_leaf = false;
	if (!q.is_leaf())
	{
		sum += firstraverse(*q.children[CHILD_LEFT], r);
		sum += firstraverse(*q.children[CHILD_RITE], r);
	}
	q.processed = true;
	if (q.is_leaf())
		return sum;


	t_kitsch_node& ii = *q.children[CHILD_LEFT];
	t_kitsch_node& jj = *q.children[CHILD_RITE];
	unsigned i = ii.get_UID();
	unsigned j = jj.get_UID();
	assert(w[i][j] == w[j][i]);
	assert(d[i][j] == d[j][i]);

	// can't recover from sibling nulls: set to impossibly large value
	if (w[i][j] <= 0.0)
	{
		q.len_to_leaf = 10.0;
		++invalid_tree;
	}
	else
		q.len_to_leaf =	d[i][j] * 0.5;

	q.weight = ii.weight + jj.weight + 2 * w[i][j];

	// gets rid of negative distances
	scrunch(q);

	sum += secondtraverse(r, q, ii, jj);


#ifndef NDEBUG
	r.CHECK_LINKS();
#endif

	return sum;
}  /* firstraverse */




//________________________________________________________________________________________

//	sumtraverse

//		traverse to finish computation of sum of squares

//________________________________________________________________________________________
float t_kitsch_tree::sumtraverse(t_kitsch_node& q)
{
	debug_func_cout(3);
	assert(!q.is_leaf());

	float sum = 0.0;

	// recurse
	if (!q.children[CHILD_LEFT]->is_leaf())
		sum += sumtraverse(*q.children[CHILD_LEFT]);
	if (!q.children[CHILD_RITE]->is_leaf())
		sum += sumtraverse(*q.children[CHILD_RITE]);


	t_kitsch_node * ii = q.children[CHILD_LEFT];
	t_kitsch_node * jj = q.children[CHILD_RITE];
	unsigned i = ii->get_UID();
	unsigned j = jj->get_UID();
	sum +=	2.0 * w[i][j] * square(d[i][j] - 2.0 * q.len_to_leaf);

	return sum;
}  /* sumtraverse */




//________________________________________________________________________________________

//	op_print_tree


//________________________________________________________________________________________
class op_check_adjacent_nulls
{
	t_kitsch_tree&		tree;
	bool				fail;
public:
	op_check_adjacent_nulls(t_kitsch_tree& tree_):
			 tree(tree_), fail(false){}
	void check(t_kitsch_node& curr_node)
	{
		vector<unsigned> UIDs;
		recurse(curr_node, UIDs);
	}
	void recurse(t_kitsch_node& curr_node, vector<unsigned>& UIDs);
	bool are_branches_linked(const vector<unsigned>& left, const vector<unsigned>& rite)
	{
		for (unsigned ii = 0; ii < left.size(); ++ii)
		{
			for (unsigned jj = 0; jj < rite.size(); ++jj)
			{
				unsigned x = left[ii];
				unsigned y = rite[jj];
				if (!tree.nulls[x][y])
					return true;
			}
		}
		return false;
	}
};

void op_check_adjacent_nulls::recurse(t_kitsch_node& curr_node, vector<unsigned>& UIDs)
{
	if (fail)
		return;

	if (curr_node.is_leaf())
	{
		UIDs.push_back(curr_node.get_UID());
		return;
	}
	vector<unsigned> branch_UIDs [CHILD__ALL];
	recurse(*curr_node.children[CHILD_LEFT], branch_UIDs[CHILD_LEFT]);
	recurse(*curr_node.children[CHILD_RITE], branch_UIDs[CHILD_RITE]);

	if (!are_branches_linked(branch_UIDs[CHILD_LEFT], branch_UIDs[CHILD_RITE]))
	{
		fail = true;
		return;
	}
	UIDs.insert(UIDs.end(), branch_UIDs[CHILD_LEFT].begin(), branch_UIDs[CHILD_LEFT].end());
	UIDs.insert(UIDs.end(), branch_UIDs[CHILD_RITE].begin(), branch_UIDs[CHILD_RITE].end());
}

//________________________________________________________________________________________

//	evaluate

//		fill in len_to_leafs and evaluate sum of squares for tree

//________________________________________________________________________________________
void t_kitsch_tree::evaluate()
{
	debug_func_cout(1);

	// cache evaluate results for this tree
	t_digest md4_hash_of_tree;
	to_md4_hash(md4_hash_of_tree);

	t_map_digest_float::iterator where = cached_evaluate_results.find(md4_hash_of_tree);
	if (where != cached_evaluate_results.end())
	{
		cnt_evaluate_cache_hits++;
		sum_sqrs = where->second;
		return;
	}

	for (unsigned i = 0; i < nodes.size(); ++i)
		nodes[i]->processed = nodes[i]->is_leaf();

/*	for (unsigned i = 0; i < nodes.size(); ++i)
		if (nodes[i]->is_leaf())
		{
			if (nodes[i]->get_parent() && nodes[i]->sibling() && nodes[i]->sibling()->is_leaf())
			{
				unsigned ii = nodes[i]->get_UID();
				unsigned jj = nodes[i]->sibling()->get_UID();
				if (w[ii][jj] == 0.0)
				{
					cerr << "!";
					sum_sqrs = -1000.0;
					cached_evaluate_results[md4_hash_of_tree] = sum_sqrs;
					return;
				}
			}
		}
*/
	invalid_tree = 0;
	float sum = firstraverse(*root, *root);
	{
		if (!root->is_leaf())
			sum += sumtraverse(*root);
		++(*pdots);
	}
	if (invalid_tree)
	{
		sum += 100000 + invalid_tree;
	}
	++cnt_trees_examined;
	sum_sqrs = -sum;
//	std::cerr << sum_sqrs << "\n";

	// save result in cache
	cached_evaluate_results[md4_hash_of_tree] = sum_sqrs;

}  /* evaluate */





//________________________________________________________________________________________

//	try_add

//		temporarily adds one fork and one leaf to the tree.
//		if the location where they are added yields greater
//		"likelihood" than other locations tested up to that
//		time, then keeps that location as there

//________________________________________________________________________________________
void t_kitsch_tree::try_add(	t_kitsch_node& p,
						   		t_kitsch_node& item,
						   		t_kitsch_node*& parent)
{
#ifndef NDEBUG
	debug_func_cout(5);
#endif


	p.left_splice_below( *parent, item);
	evaluate();
	if (sum_sqrs > best_sum_sqrs)
	{
		best_sum_sqrs = sum_sqrs;
		optimal_add_locus = &p;
	}
	parent = remove(item);



#ifndef NDEBUG
	root->CHECK_LINKS();
#endif
}  /* try_add */




//________________________________________________________________________________________

//	recursive_find_optimal_add_locus

//		traverses a binary tree, calling try_add
//		at a t_kitsch_node before calling try_add at its descendants

//________________________________________________________________________________________
void t_kitsch_tree::recursive_find_optimal_add_locus(	t_kitsch_node& p,
														t_kitsch_node& item,
														t_kitsch_node* parent)
{
#ifndef NDEBUG
	debug_func_cout(5);
#endif

	try_add(p, item, parent);
	if (!p.is_leaf())
	{
		assert(p.children[CHILD_LEFT]);
		assert(p.children[CHILD_RITE]);
		recursive_find_optimal_add_locus(*p.children[CHILD_LEFT], item, parent);
		recursive_find_optimal_add_locus(*p.children[CHILD_RITE], item, parent);
	}
}  /* recursive_find_optimal_add_locus */


//________________________________________________________________________________________

//	get_optimal_add_locus

//		gets the best place to add an item

//________________________________________________________________________________________
void t_kitsch_tree::get_optimal_add_locus(	t_kitsch_node& item,
											t_kitsch_node* parent)
{
	// set best_sum_sqrs to something ridiculous so that optimal_locus_for_add
	//		always works
	best_sum_sqrs = -1e100;
	optimal_add_locus = 0;
	recursive_find_optimal_add_locus(*root, item, parent);
	assert(optimal_add_locus);
}

//________________________________________________________________________________________

//	try_rearrange

//		evaluates one rearrangement of the tree.
//		if the new tree has greater "likelihood" than the old
//		one sets success := TRUE and keeps the new tree.
//		otherwise, restores the old tree

//________________________________________________________________________________________
void t_kitsch_tree::try_rearrange(t_kitsch_node& p)
{
#ifndef NDEBUG
	debug_func_cout(5);
#endif

	t_kitsch_node * parent		= p.get_parent();
	if (!parent)
		return;

	t_kitsch_node* grandparent	= parent->get_parent();
	if (!grandparent)
		return;

	t_kitsch_node& sibling		= *p.sibling();


	float old_sum_sqrs = sum_sqrs;

	// rearrange tree
	parent = remove(p);
	grandparent->left_splice_below(*parent, p);
	if (root->get_parent())
		root = root->get_parent();
	evaluate();

	// use that if min squares is smaller
	if (sum_sqrs > old_sum_sqrs)
	{
		best_sum_sqrs = sum_sqrs;
		sum_sqrs_lowered = true;
		return;
	}

	// else restore to previous configuration
	parent = remove(p);
	sibling.left_splice_below(*parent, p);
	if (root->get_parent())
		root = root->get_parent();
	sum_sqrs = old_sum_sqrs;

#ifndef NDEBUG
	root->CHECK_LINKS();
#endif
}  /* try_rearrange */




//________________________________________________________________________________________

//	recursive_find_optimal_rearrangement

//		traverses a binary tree, calling PROCEDURE try_rearrange
//		at a t_kitsch_node before calling try_rearrange at its descendants */

//________________________________________________________________________________________
void t_kitsch_tree::recursive_find_optimal_rearrangement(t_kitsch_node& item)
{
#ifndef NDEBUG
	debug_func_cout(5);
#endif

	try_rearrange(item);
	if (!item.is_leaf())
	{
#ifndef NDEBUG
		debug_cout(5, item.get_UID());
		debug_cout(5, item.children[CHILD_LEFT]->get_UID() << "left");
		debug_cout(5, item.children[CHILD_RITE]->get_UID() << "right");
#endif
		recursive_find_optimal_rearrangement(*item.children[CHILD_LEFT]);
		recursive_find_optimal_rearrangement(*item.children[CHILD_RITE]);
	}
}  /* recursive_find_optimal_rearrangement */




//________________________________________________________________________________________

//	rearrange

//		traverses the tree (preorder), finding any local
//		rearrangement which decreases the number of steps.
//		if traversal succeeds in increasing the tree's
//		"likelihood", PROCEDURE rearrange runs traversal again

//________________________________________________________________________________________
void t_kitsch_tree::rearrange()
{
#ifndef NDEBUG
	debug_func_cout(5);
#endif


	sum_sqrs_lowered = false;
	do
	{
		sum_sqrs_lowered = false;
		recursive_find_optimal_rearrangement(*root);
	} while (sum_sqrs_lowered);
}  /* rearrange */





t_kitsch_tree& 	t_kitsch_tree::operator = (const t_kitsch_tree& other)
{
	// these are scratch variables
	optimal_add_locus = 0;
	sum_sqrs_lowered = false;

	assert(cnt_leaves		== other.cnt_leaves);
	assert(&names			== &other.names);
	assert(&w				== &other.w);
	assert(&d				== &other.d);
	assert(pdots			== other.pdots);
	assert(power			== other.power);
//	names					= other.names;
//	cnt_leaves	 				= other.cnt_leaves;
//	cnt_nodes				= other.cnt_nodes;
	cnt_trees_examined		= 0;//other.cnt_trees_examined;
	cnt_non_null_data		= other.cnt_non_null_data;
	shared_mem				= other.shared_mem;
	sum_sqrs                = other.sum_sqrs;
	best_sum_sqrs           = other.best_sum_sqrs;

	nodes.assign(cnt_nodes, 0);
	for (unsigned i = 0; i < cnt_nodes; ++i)
		nodes[i] = &shared_mem[i];

	for (unsigned i = 0; i < cnt_nodes; ++i)
	{
		if (other.nodes[i]->get_parent())
		{
			// restore link to parent
			nodes[i]->parent = nodes[other.nodes[i]->get_parent()->treepos];
			nodes[i]->get_parent()->children[nodes[i]->child_num] = nodes[i];
		}
		else nodes[i]->parent	= 0;
	}
	if (other.root)
		root = nodes[other.root->treepos];
	else
		root = 0;
	return *this;
}


//________________________________________________________________________________________

//	maketree


//		constructs a binary tree from the pointers in treenode.
//		adds each t_kitsch_node at location which yields highest "likelihood"
//		then rearranges the tree for greatest "likelihood"

//________________________________________________________________________________________
void t_kitsch_tree::maketree(const vector<unsigned>& order_of_entry)
{
#ifndef NDEBUG
	debug_func_cout(2);
#endif

	for (unsigned i = 0; i < cnt_nodes; ++i)
		nodes[i]->clear();

	// first two leaves
	root = nodes[order_of_entry[0]];
	nodes[order_of_entry[0]]->left_splice_below(*nodes[cnt_leaves], *nodes[order_of_entry[1]]);

	// add all other leaves
	for (unsigned i = 2; i < cnt_leaves; ++i)
	{
		// get next leaf and the next unused internal node
		t_kitsch_node& item			= *nodes[order_of_entry[i]];
		t_kitsch_node* new_branch	= nodes[cnt_leaves - 1 + i];


		// add the node in the best place
        get_optimal_add_locus(item, new_branch);
		optimal_add_locus->left_splice_below(*new_branch, item);

		// after adding in the optimal locus, the sum_sqrs should also be optimal:
		sum_sqrs = best_sum_sqrs;

		// rearrange everything locally to see if can get a better optimal
		rearrange();

		// evaluate without bumping up the number of trees examined
		evaluate();
		--cnt_trees_examined;
	}



	// take each node out and add back in at the optimal position
	// keep doing until this doesn't improve the sum of squares any more
	float orig_sum_sqrs;
	do
	{
		orig_sum_sqrs = best_sum_sqrs;
		t_kitsch_node * new_branch = 0;
		for (unsigned j = 0; j < cnt_nodes; ++j)
		{
			t_kitsch_node& item = *nodes[j];
			if (&item != root)
			{
				new_branch = remove(item);
				get_optimal_add_locus(item, new_branch);
				optimal_add_locus->left_splice_below(*new_branch, item);
			}
		}
	} while (best_sum_sqrs > orig_sum_sqrs);

	debug_cout(4, "best_sum_sqrs " << best_sum_sqrs);
}  /* maketree */


void t_kitsch_tree::maketree(const bcpp_tree::t_tree& orig_tree)
{
#ifndef NDEBUG
	debug_func_cout(2);
#endif

	if (orig_tree.get_node_count() != cnt_nodes)
		throw std::logic_error("Different number of nodes between Kitsch tree (" +
								ultoa(cnt_nodes) + ") and original tree (" +
								orig_tree.get_name() + " " +
								ultoa(orig_tree.get_node_count()) + ")");
	for (unsigned i = 0; i < cnt_nodes; ++i)
		nodes[i]->clear();


	bcpp_tree::t_tree normalized_tree;
	orig_tree.get_normalized_tree(normalized_tree);
	root = nodes[normalized_tree.get_root()->get_UID()];
	for (unsigned i = 0; i < normalized_tree.get_node_count(); ++i)
	{
		const bcpp_tree::t_node& orig_node = normalized_tree.get_node(i);
		if (orig_node.get_parent() != 0)
		{
			nodes[i]->parent = nodes[orig_node.get_parent()->get_UID()];
		}
		if (!orig_node.is_leaf())
		{
			assert(!nodes[i]->leaf);
			nodes[i]->children[0] = nodes[orig_node.get_child(0)->get_UID()];
			nodes[i]->children[1] = nodes[orig_node.get_child(1)->get_UID()];
			nodes[i]->children[0]->child_num = kitsch_tree::CHILD_LEFT;
			nodes[i]->children[1]->child_num = kitsch_tree::CHILD_RITE;
		}
	}
	CHECK_LINKS();
	// evaluate without bumping up the number of trees examined
	evaluate();
	--cnt_trees_examined;

	// rearrange everything locally to see if can get a better optimal
	rearrange();

	// evaluate without bumping up the number of trees examined
	evaluate();
	--cnt_trees_examined;



	// take each node out and add back in at the optimal position
	// keep doing until this doesn't improve the sum of squares any more
	float orig_sum_sqrs;
	do
	{
		orig_sum_sqrs = best_sum_sqrs;
		t_kitsch_node * new_branch = 0;
		for (unsigned j = 0; j < cnt_nodes; ++j)
		{
			t_kitsch_node& item = *nodes[j];
			if (&item != root)
			{
				new_branch = remove(item);
				get_optimal_add_locus(item, new_branch);
				optimal_add_locus->left_splice_below(*new_branch, item);
			}
		}
	} while (best_sum_sqrs > orig_sum_sqrs);

	debug_cout(4, "best_sum_sqrs " << best_sum_sqrs);
}  /* maketree */












//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	treeout

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
template<typename TREEMAKER> void kitsch_treeout_recurse(	TREEMAKER& maker,
															t_kitsch_tree& tree,
															t_kitsch_node *p,
															const vector<string>& names
															)
{
#ifndef NDEBUG
	debug_func_cout(2);
#endif

	assert(p);
	if (p->is_leaf())
	{
		maker.add_leaf(names[p->get_UID()], p->get_UID());
	}
	else
	{
		maker.add_internal_node();
		kitsch_treeout(maker, tree, p->children[CHILD_LEFT], names);
		maker.add_sibling();
		kitsch_treeout(maker, tree, p->children[CHILD_RITE], names);
		maker.reparent_to_grandparent();
	}
	if (p == tree.root)
	{
		maker.finish();
	}
	else
	{
		maker.set_branch_len(p->get_parent()->len_to_leaf - p->len_to_leaf);
	}
}

template<typename TREEMAKER> void kitsch_treeout(	TREEMAKER& maker,
													t_kitsch_tree& tree,
													t_kitsch_node *p,
													const vector<string>& names
													)
{
#ifndef NDEBUG
	debug_func_cout(2);
#endif

	assert(p);
	if (p->is_leaf())
	{
		maker.add_leaf(names[p->get_UID()], p->get_UID());
	}
	else
	{
		maker.add_internal_node();
		kitsch_treeout(maker, tree, p->children[CHILD_LEFT], names);
		maker.add_sibling();
		kitsch_treeout(maker, tree, p->children[CHILD_RITE], names);
		maker.reparent_to_grandparent();
	}
	if (p == tree.root)
	{
		maker.finish();
	}
	else
	{
		maker.set_branch_len(p->get_parent()->len_to_leaf - p->len_to_leaf);
	}
}

unsigned t_kitsch_tree::recurse_reorder(unsigned parent_UID) const
{
	debug_func_cout(1);
	assert(parent_UID < normalisation_nodes.size());
	t_reorder_node& node = normalisation_nodes[parent_UID];
	if (node.is_leaf)
		return node.leaf_UID;

	node.leaf_UID = recurse_reorder(node.left);
	unsigned rite_UID = recurse_reorder(node.rite);
	if (rite_UID < node.leaf_UID)
	{
		std::swap(node.left,  node.rite);
		node.leaf_UID = rite_UID;
	}
	return node.leaf_UID;
}
/*

struct t_state{
};
unsigned t_kitsch_tree::reorder(unsigned parent_UID) const
{
	debug_func_cout(1);
	assert(parent_UID < normalisation_nodes.size());
	t_reorder_node& node = normalisation_nodes[parent_UID];
	if (node.is_leaf)
		return node.leaf_UID;
	deque<t_state> call_stack;
	while (call_stack.size())
	{
		unsigned left_UID = recurse_reorder(node.left);
		unsigned rite_UID = recurse_reorder(node.rite);
		if (rite_UID < left_UID)
		{
			std::swap(node.left,  node.rite);
			node.leaf_UID = rite_UID;
		}
		else
			node.leaf_UID = rite_UID;
		return node.leaf_UID;
	}
}

*/



t_digest& t_kitsch_tree::to_md4_hash(t_digest& hv) const
{
	debug_func_cout(1);
	normalisation_nodes.assign(shared_mem.begin(),  shared_mem.end());
	assert(root);
	unsigned root_UID = root->UID;
	recurse_reorder(root_UID);
	vector<unsigned> left_rite(normalisation_nodes.size() * 2, 0);
	for (unsigned i = 0; i < normalisation_nodes.size(); ++i)
	{
		left_rite[i * 2] = normalisation_nodes[i].left;
		left_rite[i * 2 + 1] = normalisation_nodes[i].rite;
	}

	debug_cout(1, "now to md4 hash");
	return md4_hash(reinterpret_cast<unsigned char*>(&left_rite[0]),
					left_rite.size() * sizeof(unsigned), hv);
}




string t_kitsch_tree::to_str()
{
	newhampshire_tree_maker make_newhampshire_tree(true);
	kitsch_treeout(make_newhampshire_tree, *this, root, names);
	return make_newhampshire_tree.str;
}













//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	use kitsch to make a paralogy tree

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
void kitsch_make_tree(	bcpp_tree::t_tree&				tree,
						const std::vector<std::string>&	names,
						const bcpp_matrix::t_dense_matrix_f&	matrix,
						unsigned						njumble,
						float							power,
						unsigned&						best_jumble,
						float&							best_sum_sqrs,
						float&							avg_percent_stddev,
						unsigned&						cnt_trees_examined,
						t_progress_indicator&			dots)
{
#ifndef NDEBUG
	debug_func_cout(1);
#endif
	best_jumble = 1;
	unsigned cnt_leaves = matrix.size();
	unsigned cnt_nodes = cnt_leaves * 2 - 1;

	bcpp_matrix::t_dense_matrix_f  		weights(cnt_nodes, 0.0);
	bcpp_matrix::t_dense_matrix_f		matrix_copy(cnt_nodes, 0.0);
	matrix_copy.copy_from(matrix);

	t_kitsch_tree curr_tree(matrix_copy, weights, names, power, dots);
	// initialize distances and weights

	curr_tree.prepare_data();
	t_kitsch_tree best_tree(matrix_copy, weights, names, power, dots);

	// which order leaves are added
	vector<unsigned> order_of_entry;
	for (unsigned i = 0; i < names.size(); ++i)
		order_of_entry.push_back(i);


	//
	// make initial upgma tree
    //
    bcpp_matrix::t_sparse_matrix_f  sparse_matrix(matrix, 10000.0);
	bcpp_tree::t_tree upgma_tree;
	upgma_tree::upgma_make_tree(upgma_tree, sparse_matrix,false);
	curr_tree.maketree(upgma_tree);
	if (curr_tree.best_sum_sqrs > best_tree.best_sum_sqrs)
	{
		best_tree = curr_tree;
		best_jumble	= 0;
	}
	cnt_trees_examined += curr_tree.cnt_trees_examined;


	for (unsigned i = 1; i < njumble + 1; ++i)
	{
		// jumble
		std::reverse(order_of_entry.begin(), order_of_entry.end());
		for (unsigned j = 0; j < order_of_entry.size(); ++j)
			std::swap(order_of_entry[j], order_of_entry[rand() % order_of_entry.size()]);

		curr_tree.maketree(order_of_entry);
		if (curr_tree.best_sum_sqrs > best_tree.best_sum_sqrs)
		{
			best_tree = curr_tree;
			best_jumble	= i + 1;
		}
		cnt_trees_examined += curr_tree.cnt_trees_examined;
	#ifndef NDEBUG
			best_tree.CHECK_LINKS();
			curr_tree.CHECK_LINKS();
	#endif

	}

	best_tree.evaluate();

	best_sum_sqrs		= -best_tree.sum_sqrs;
	avg_percent_stddev	= 100.0 * sqrt(best_sum_sqrs /
									(best_tree.cnt_non_null_data - 1.0) / 2.0);

	// kitsch tree => paralogy tree
	bcpp_tree_maker make_bcpp_tree(tree);
	kitsch_treeout(make_bcpp_tree, best_tree, best_tree.root, names);
	dots.finish();
}

//________________________________________________________________________________________

//		kitsch_make_tree

//			wrapper

//________________________________________________________________________________________
void kitsch_make_tree(	bcpp_tree::t_tree&			tree,
						const std::vector<std::string>&	names,
						const bcpp_matrix::t_dense_matrix_f&	matrix,
						unsigned						njumble,
						float							power,
						unsigned&						best_jumble,
						float&							sum_sqrs,
						float&							avg_percent_stddev,
						unsigned&						cnt_trees_examined)
{
#ifndef NDEBUG
	debug_func_cout(1);
#endif

	ostringstream ignore_strm;
	t_progress_indicator dots(ignore_strm, 10000);
	kitsch_make_tree(tree, names, matrix, njumble, power, best_jumble,
					sum_sqrs, avg_percent_stddev,
						cnt_trees_examined, dots);
}









//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	make a kitsch tree from a paralogy tree

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//________________________________________________________________________________________

//	do_create_kitsch_tree

//________________________________________________________________________________________
struct do_create_kitsch_tree
{
	do_create_kitsch_tree(  kitsch_tree::t_kitsch_tree&				new_tree_,
							bcpp_tree::t_tree					tree_		)

						:	new_tree(new_tree_),
							cnt_leaves(new_tree_.cnt_leaves), current_node(0),
							next_internal_node(new_tree.cnt_leaves),
							next_seq_node(0)
	{
		new_tree.root = new_tree.nodes[next_internal_node];
		bcpp_tree_ops::do_calculate_branch_lengths branch_lengths(tree_);
		lengths_to_leaves.swap(branch_lengths.lengths_to_leaves);
	}

	std::vector<float>	 					lengths_to_leaves;
	kitsch_tree::t_kitsch_tree&				new_tree;
	unsigned								cnt_leaves;
	kitsch_tree::t_kitsch_node*				current_node;
	unsigned								next_internal_node;
	unsigned								next_seq_node;
	void									operator()(bcpp_tree::t_tree::t_node_type* root)
												{recurse(root, CHILD_LEFT);}
	void									recurse(bcpp_tree::t_tree::t_node_type* node, eCHILD_NUM child_num);
};


//________________________________________________________________________________________

//	do_create_kitsch_tree

//________________________________________________________________________________________
void do_create_kitsch_tree::recurse(bcpp_tree::t_tree::t_node_type* node, eCHILD_NUM child_num)
{
	assert(node);

	// create leaf
	if (node->is_leaf())
	{
		assert(new_tree.nodes[next_seq_node]->is_leaf() == true);
		assert(next_seq_node < new_tree.cnt_leaves);
		new_tree.nodes[next_seq_node]->len_to_leaf = lengths_to_leaves[node->get_UID()];
		new_tree.nodes[next_seq_node]->UID = node->get_leaf_UID();
		new_tree.nodes[next_seq_node++]->set_parent(current_node, child_num);
		return;
	}


	// add_internal_node
	assert(new_tree.nodes[next_internal_node]->is_leaf() == false);
	assert(next_internal_node < new_tree.cnt_nodes);

	new_tree.nodes[next_internal_node]->len_to_leaf = lengths_to_leaves[node->get_UID()];
	new_tree.nodes[next_internal_node]->set_parent(current_node, child_num);
	current_node = new_tree.nodes[next_internal_node];
	++next_internal_node;

	recurse(node->get_child(0), CHILD_LEFT);
	recurse(node->get_child(1), CHILD_RITE);

	current_node = current_node->get_parent();
}















//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	Calculate sum of squares for a paralogy tree

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
void optimize_bcpp_tree(	bcpp_tree::t_tree&				tree,
							const vector<string>&			names,
							const bcpp_matrix::t_dense_matrix_f&	matrix,
							unsigned						njumble,
							float							power,
							unsigned						max_at_a_time,
							unsigned&						best_jumble,
							float&							best_sum_sqrs,
							float&							avg_percent_stddev,
							unsigned&						cnt_trees_examined,
							t_progress_indicator&			dots)
{
	best_jumble = 0;

	// if tree too small just kitsch it all
	if (max_at_a_time >= tree.get_leaf_count())
	{
		tree.reset();
		kitsch_make_tree(tree, names, matrix, njumble, power, best_jumble,
							best_sum_sqrs, avg_percent_stddev,
							cnt_trees_examined, dots);
		return;
	}
//	else
//	{
//		cerr << max_at_a_time << ":" <<  tree.get_leaf_count() << std::endl;
//	}

	// allocate memory for distances and weights
	unsigned		cnt_nodes = matrix.size() * 2 - 1;
	bcpp_matrix::t_dense_matrix_f  		weights(cnt_nodes, 0.0);
	bcpp_matrix::t_dense_matrix_f		matrix_copy(cnt_nodes, 0.0);
	matrix_copy.copy_from(matrix);


	// make initial kitsch tree from paralogy tree
	t_kitsch_tree optimized_tree(matrix_copy, weights, names, power, dots);
	do_create_kitsch_tree tree_maker(optimized_tree, tree);
	optimized_tree.prepare_data();
	tree.for_each(tree_maker);
	optimized_tree.evaluate();

	// optimize sub trees using kitsch
	op_optimize_sub_trees(optimized_tree, njumble, max_at_a_time, cnt_trees_examined).
							start_at(*optimized_tree.root);
	best_sum_sqrs		= -optimized_tree.sum_sqrs;
	avg_percent_stddev	= 100.0 * sqrt(best_sum_sqrs /
										(optimized_tree.cnt_non_null_data - 1.0) / 2.0);


	// kitsch tree => paralogy tree
	tree.reset();
	bcpp_tree_maker make_bcpp_tree(tree);
	kitsch_treeout(make_bcpp_tree, optimized_tree, optimized_tree.root, names);

}





} // namespace

